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NUMERICAL SIMULATION OF PARABOLIC 
MOVING AND GROWING INTERFACE PROBLEMS 
USING SMALL MESH DEFORMATION 


ULRICH LANCER AND HUIDONG YANG 

Abstract. In this work, we develop a cutting method for solv¬ 
ing problems with moving and growing interfaces in 3D. This new 
method is able to resolve large displacement or deformation of im¬ 
mersed objects by combining the Arbitrary Lagrangian-Eulerian 
method with only small local mesh deformation defined on the 
reference domain, that is decomposed into the macro-elements. 
The linear system of algebraic equations arising after the tempo¬ 
ral and spatial discretizations of a model parabolic interface heat- 
conduction-like problem with vector-valued functions is solved by 
either an all-at-once or a segregated algebraic multigrid method. 


1. Introduction 

The conventional Arbitrary Lagrangian-Eulerian (ALE) method (see, 
e.g., Ha El) works well for small deformation in many applications. 
For large deformation problem, the ALE method may fail due to the 
deteriorated mesh quality. Some improved ALE methods have been 
studied, e.g., a method based on the biharmonic extension in |2T] . 
More work based on the so called fixed-mesh ALE approach has been 
studied, e.g., in p U- The parametric finite element method [8], 
the immersed-interface finite element method na and the immersed 
boundary method [2] may also be applied in this context. An enhanced 
ALE method combined with the fixed-grid and extended finite element 
method (XEEM) was studied in [Uj. Another promising approach is 
to use the space-time method, that is more flexible to handle moving 
interface problems; see, e.g., dZlIIH]. 

In this work, we propose an interface capturing method by pre- 
computing the intersection of the moving object immersed in the un¬ 
derlying reference tetrahedral elements in three dimension (3D). Com¬ 
bined with the ALE method on such reference elements, we are able 
to deal with the moving or growing interface problems with large dis¬ 
placement or deformation. In a similar manner as already investigated 
in the earlier work [HI EH E2] , the piece-wise linear finite element basis 
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functions are constructed on each macro-element [19], that is decom¬ 
posed into four pure tetrahedral elements and one octahedral element. 
In addition, the method offers a nice opportunity to keep capturing 
the interface without introducing extra degrees of freedom. To test the 
robustness of the method, we consider a model heat-conduction-like 
problem with vector-valued functions. Such a model can be used to 
handle the mesh movement in the fluid-structure interaction simula¬ 
tion, see, e.g., [TH]. The construction of robust solution methods for 
solving the arising finite element equations requires additional effort. 
For this, we use both the all-at-once and the segregated methods, that 
employ an algebraic multigrid (AMG) method [TTl fTT] . 

The remainder of the paper is organized as follows: In Section we 
set up the model parabolic interface problem. Section deals with the 
temporal and spatial discretization of the model interface problem. In 
Section |4| we discuss the all-at-once and the segregated methods for 
solving the linear system of equations arising from the temporal and 
spatial discretization. We present numerical results of two proposed 
interface moving problems in Section Finally, some conclusions are 
drawn in Section [6l 

2. A MODEL INTERLACE PROBLEM 

2.1. Geometrical configurations. We consider a simply connected, 

bounded, polyhedral Lipschitz domain G C which includes an im¬ 
mersed time-dependent, sufficiently smooth sub-domain fl\ = fliit) C 
fl, where t £ I denotes the time with I = (0, T] being the time inter¬ 
val. The remaining sub-domain is = r2\f25- By F* := dfl\ we 

denote the interface. The boundaries of are denoted by Fd and Fat 
such that dVL = f £> U Fjv and F^) n F^r = 0, where proper Dirichlet and 
Neumann boundary conditions are prescribed, respectively. We use n 
to denote the outward unit normal vector on the boundary dfl, ni and 
n 2 , the outward unit normal vectors on F* with respect to Gi and G 2 , 
respectively. We refer to Fig. [T]for an illustration of a sub-domain im¬ 
mersed in the big domain. We consider two interface problems in this 
work. In the first problem, the sub-domain keeps the shape and 
moves with a constant velocity n G i.e., a rigid body motion; see 
the left plot in Fig. In the second problem, the sub-domain grows 
with a constant velocity v along the line connecting the mass center 
of the sub-domain and any point pt on the boundary see the right 
plot in Fig. [1} 

2.2. The model problem with a fixed interface F. We start to 
formulate the problem in the fixed sub-domains and O 2 with proper 
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Figure 1. An illustration of two sub-domains for the 
interface problem: Rigid body motion (left) and growing 
sub-domain (right). 


interface conditions on the fixed interface F = fl d^ 2 - We aim to 
find the solution n : fl i-?- for alH G I, such that 


( 1 ) 


dfU — V • (aVn) = 0 in U ^ 2 , 

Ui = U 2 on r, 

dui du2 ^ 
ai— - h a 2 -— on i 


drii 


‘dn2 


with the initial condition n = 0 at t = 0, and the boundary conditions 
U 2 — gn onT]:, and Q- 2 |^ = (/jv on F^r at t > 0. Here a = Oi G M"*" in fli, 
a = 02 G in fl 2 , ai a 2 , are two different material coefficients. The 
analysis of such an interface problem with the scalar-valued function 
has been studied, e.g., in 0 HSl H- In this work, we consider the 
model problem with the vector-valued function, that can be used to 
model the mesh movement in the fluid-structure interaction simulation 
in our future work. 


2.3. The model problem with a unfixed interface Fh For the 

interface problem with unfixed interface F*, the time derivative dtu in 
Q is not well-defined since the computational domain is moving. One 
of the classical approaches is to use the ALE method iniiz], in which 
we introduce a displacement defined on the reference domain Qr: 


d{x, t) : flji X / i-A 
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for all X G Vtji and t ^ that tracks the monition of the computational 
domain Q. The ALE mapping ^ for all t G /, where 

= Cl\U ^ 2 ^ is defined as 

= A{x^ t) \= X + d{x^ t) 

for all X G VLr and t ^ I. In our model problem, we shall interprete d 
as the finite element mesh movement, which defines the change of the 
computational sub-domains and is explicitly precomputed. The ALE 
time derivate of the function r : i-A is defined as 


dtu\Ai := dtu{A\x,t),t) 

for all X G flR and t £ I. By the chain rule, we obtain 


dtu = dtu\jxt — w • Vn 

with w = dfA^ o A^~^. Then we have the following model problem 
under the ALE framework: Find the solution u : i—>■ for all 

t £ I, such that 


( 2 ) 


~ w ■ Vu — V • (aVu) = 0 in n* U 

ui = U 2 on C, 

dui du2 
Cl\~ - CL2~ - DU i 


drii 


‘dm 


with the initial conditions u{x, 0) = 0, w{x^ 0) = 0 for all a: G U Q®) 
and the boundary conditions u = gr, on Tr, and a|^ = pjv on Ttv at 
t > 0. Here a = ai G M"*" in a = a 2 £ M'*' in 0^, Ui ^ a 2 , are two 
different material coefficients in two moving domains, respectively. 


2.4. A combination of the ALE and macro-element method. In 

the classical ALE method, we use the interface tracking method, where 
the mesh grids on the interface are following the object movement. The 
mesh movement inside the computational domains is computed by an 
arbitrary extension into the domain, e.g., a simple harmonic exten¬ 
sion. The main drawback of this method is the restriction to small 
deformations. In case of large deformation or displacement, the mesh 
quality may deteriorate rapidly. To overcome this difficulty, we develop 
an interface capturing method, that is a combination of the ALE and 
macro-element method [T9l Eg. According to the cutting cases, the 
underlying reference domain is decomposed into macro-elements: four 
triangles in each macro-element in 2D and four tetrahedra plus one 
octahedron in 3D, see Fig. [^for an illustration of such decomposed ref¬ 
erence domain into structured grids in 2D. The velocity w : flji 
of the mesh movement is constructed locally in each sub-element of the 
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macro-element by an interpolation. The same applies to the displace¬ 
ment d : Qr I—>■ of the mesh movement, with respect to the reference 
configuration Qr. The local velocity and displacement are related by 
w = dtd. We comment that, for cells that are completely untouched 
with the moving interface (far away from the moving object), the ve¬ 
locity tc = 0 and the ALE mapping is an identity. In this case, the 
equation ([^ is reduced to the one under the usual Eulerian framework. 





Eigure 2. An illustration of a reference domain Q,r 
decomposed into macro-elements: macro-element edges 
(thick solid lines), introduced new sub-element edges 
(thin dashed lines) and the interface (thick blue dashed 
line). 


3. Temporal and spatial discretization 

3.1. Temporal discretization. Let the time interval I be divided 
into N equidistant small time intervals At, i.e.. At = T/N. Let 
= nAt be the time at level n. By the notation /” = f(x,t^), 
we denote the function defined at the time t” and in the corresponding 
domain. We employ first-order implicit Euler scheme to discretize the 
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Figure 3. An illustration of the local nodes movement 
with restriction to each macro-element: hxed macro¬ 
element nodes (brown dots), reconstructed moving in¬ 
terface and locally adapted triangle mesh at t = with 
nodes from the intersection (red lines), reconstructed 
moving interface and locally adapted triangle mesh at 
t = with nodes from the intersection (blue lines), the 
moving direction of the intersection nodes within each 
macro-element at the interface, cyan arrows (none of the 
intersection nodes is the edge middle point), magenta 
arrows (one of the nodes is the edge middle point). 


time derivative: For all n > 1 and given = uq, we have 


u — u 


n—1 


At 


— w 


Vn” - V • (aVn”) = 0 in U flf, 


( 3 ) 


< = «" on F*", 
dvtt duo _|ri 
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where o o te” = o In our 

model problem, the displacement d is dehned on the reference domain 
Qr and is explicitly evaluated by the intersection of the moving object 
with the underlying tetrahedral mesh. In Fig. the cyan and magenta 
arrows indicate the movement of the intersection points with respect to 
the underlying reference macro-element mesh. Assume the underlying 
mesh size consisting of the macro-element is h, then the mesh movement 
velocity is controlled by |?n| ~ When we choose sufficiently small 
mesh size, it will only introduce a very small convection term in the 
model problem, that is in general less problematic to perform standard 
finite element discretization and to solve the arising linear system of 
equations. 

3.2. Spatial discretization. The weak formulation arises from 0 by 
integration by parts and reads as follows: Find the solution u'^ &Vg = 
{qd + bo} with Vo = = {n G H^{QY\v = 0 on F^)} such that, 

for all n G Vo, we have 

( 4 ) 

_ ^n-1 \ 

^ - ,vj ^ + (aVu”, = {gN,v)rj,, 

where the continuity condition for the solution on F* has been explicitly 
enforced by using one identical in the domain fl*”, and the surface 
traction balance condition is implicitly included in the week form by 
integration by parts. 

We use a finite element method for the spatial discretization. This 
method relies on the piecewise linear basis functions constructed on 
the underlying hybrid mesh consisting of tetrahedral and octahedral 
elements. Such mixed elements are obtained by decomposing each 
macro-element (a big tetrahedra) into four tetrahedral elements and 
one octahedral element; see Fig. for an illustration of such a typ¬ 
ical macro-element. Each tetrahedral macro-element has four fixed 
nodes with local node numbering 0 — 3 (brown dots in Fig. and 
six nodes 4 — 9 on edges (cyan dots in Fig. that are given by the 
edge middle points or the intersection points between the edge and 
the moving object. Each macro-element is decomposed into five sub¬ 
elements: four tetrahedron with the local node numbering {0,4, 6,7}, 
(4,1,5,6}, (6, 5, 2, 9}, {7,8,9,3} and one octahedron with the local 
node numbering {6,4,5,9, 7, 8}. This gives a very limited intersection 
patterns. In addition every macro-element has very similar structure to 
each other, that is easy to templatize on the computer implementation. 
By this means, we are able to reconstruct the triangle surface mesh of 
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the immersed object; see Fig. for an illustration of a sequence of such 
surface meshes. The hybrid mesh, consisting of different element types, 
has also been used recently in the cardiac electrophysiology simulation 
[19] and in the fluid-structure interaction simulation [21 Eg. 



Figure 4. An illustration of a macro-element with five 
small sub-elements. 



Figure 5. A sequence of reconstructed surface meshes 
of the immersed growing objects. 

To be more precise, the finite element basis functions on the four 
tetrahedra in each macro-element is constructed as the standard hat 
function in 3D. On the remaining octahedron, we first add an auxiliary 
point 6 near or at the mass center. The octahedron will be sub-divided 
into 8 tetrahedra; see Fig.j^for an illustration. We then construct stan¬ 
dard hat functions on each tetrahedron. The extra degree of freedom 
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at the node 6 will be eliminated by the averaging of the values at nodes 
0 — 5; see more details in dSlEl]. By this means, we do not introduce 
new degrees of freedom. The number of total degrees of freedom is the 
number of nodes plus edges in the original mesh consisting of pure big 
tetrahedral macro-elements. 


5 



Figure 6. Splitting of an octahedron into 8 tetrahedra 
(0,1,2,6}, {0,2,3,6}, (0,3,4,6}, (0,1,4,6}, (5,1,4,6}, 
(5,1,2,6}, {5,2,3,6}, (5,3,4,6}: Original edges (thick 
lines), added edges (thin lines), original nodes (0 — 5}, 
added node {6}. 


4. Solution methods eor the linear system oe equations 

4.1. An all-at-once method. After using finite element discretiza¬ 
tion, at each time step, we obtain the following linear system of equa¬ 
tions: 


' Avv 

Ave 


Uy 


' fv' 

Aev 

Aee 


Ue 


. . 


We solve the linear system of equations by the AMG preconditioned 
conjugate gradient (PCG) method (see, e.g., [H]) and the AMG pre¬ 
conditioned GMRES method (see [20]). We mention here that, due to 
the small convection term, we found out that even the PGG method 
works well for solving such a non-symmetrically perturbed symmetric 
linear system of equations. For convenience of the solution procedure, 
the linear system has been ordered with firstly the degrees of freedom 
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on the original tetrahedral nodes uy, and then of the edges ue, where 
the subscripts V and E are associated with the nodes and edges. Such 
reordering has been used in the AMG method for high-order finite el¬ 
ement discretized equations |23l [I5] . The stiffness matrices Ayv and 
Aee arise from the finite element assembly of the basis functions asso¬ 
ciated with the original macro-element nodes and edges, respectively, 
Agy and AyE are coming from the coupling. To solve such a linear 
system of equations, we use a special AMG method [Hj, that is based 
on the matrix graph connectivity. Similar idea was also developed in 
j3]. In our numerical simulation, such solution methods give us quite 
satisfactory results. We observe a quite robust behavior of the AMG 
preconditioner with respect to moving interface in each time step. 

4.2. A segregated method. By a close look at the matrix structure 
in (§, we have observed that Ayy is a block-diagonal matrix. This is 
due to the fact that the degrees of freedom associated with the original 
macro-element nodes are completely decoupled. In Fig. we demon¬ 
strate a sparsity pattern of the system matrix K, where it is easy to 
see the block-diagonal structure of Ayy. 



Figure 7. Sparsity pattern of the system matrix K. 


We now 

(S) 


perform a LU factorization of the system matrix K in 

I 


K = LU = 


’ Ayy 

0 ■ 


Aev 

s 



0 


Ayy^VE 

I 


(§: 
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where S denotes the Schur complement S = Aee — AEvAyyAvE- 
Since Ayy can be constructed very easily, the Schur complement S can 
also be constructed exactly. A simple blockwise forward and backward 
substitution gives rise to the solution of the linear system. The main 
cost is to solve the Schur complement equation 

(7) Sxee = bE 

for := /e — Aev Ayy fv ■ This is realized by applying the AMG 
preconditioned CG method m- 

5. Numerical results 

5.1. The numerical result for the model problem with an im¬ 
mersed moving sphere. In the hrst example, we consider a sphere 
with fixed radius 0.12 and the initial center at (0.125,0.125,0.125) im¬ 
mersed in a unit cube; see Fig. for an illustration. The cube is 



Figure 8. Gutting plane (left), constructed moving 
sphere surfaces at the time t — 0 (middle) and t = 0.5625 
(right). 

decomposed into macro-elements with 35937 nodes and 196608 tetra- 
hedra. The sub-divided hybrid mesh consists of 274625 nodes and 
786432 tetrahedra and 196608 octahedra. The total number of degrees 
of freedom is 823875. See Fig. [^for an illustration. 

The sphere is moving along the line with the starting point (0,0,0) 
and the ending point (1,1,1), and the moving speed is v = (1,1,1)^. 
The constructed sphere surface is shown in the middle and right plots 
of Fig. 1^ On the bottom of the cube, we set the Dirichlet boundary 
condition u = (0,0,0)^, on the top, u = (1, 0, 0)^. For the rest of the 
boundaries, we use the homogeneous Neumann boundary condition. 
The time stepsize is At — 0.0625 and the number of time steps is 9, i.e., 
the ending time is T = 0.5625. The material coefficient inside G* is Ui = 
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Figure 9. Original pure tetrahedral mesh (left), the 
sub-divided mesh (right). 


l.Oe -h 06 and inside Vt\ is ai = 1.0. The simulation results at different 
time on the cutting plane (see the left plot in Fig. are shown in 
Fig.[T^ The relative residual error is set to lOe—09 as stopping criteria. 
The iteration numbers and the computational CPU time (in second) 
of the AMG preconditioned CG and GMRES methods in the all-at- 
once method, and the iteration numbers of the AMG preconditioned 
CG for the Schur complement equation and the computational GPU 
time in the segregated method, are shown in Fig. [m We observe that, 
in terms of iteration numbers, the GMRES method shows the best 
performance, then the CG method, and last the segregated method. 
However, regarding CPU time, we see that, the CG method shows its 
best performance, then the segregated method, and last the GMRES. 


5.2. The numerical result for the model problem with au im¬ 
mersed growiug sphere. In the second example, we consider a sphere 
with an initial radius 0.08 and the initial center at (0.5,0.5, 0.5) im¬ 
mersed in a unit cube; see Fig. for an illustration. We use the same 
finite element mesh as in the first example. The sphere is growing along 
the radius direction and the growing speed is v = n, where n denotes 
the outward unit normal vector in the radius direction. The surfaces of 
the growing sphere at time t = 0 and t = 0.45 are constructed as shown 
in the middle and right plots of Fig. respectively. On the bottom 
of the cube, we set the Dirichlet boundary conditions u = (0, 0,0)^, on 
the top, u = (1,0, 0)^. For the rest of the boundaries, we use the homo¬ 
geneous Neumann boundary condition. The time stepsize is At = 0.05 
and the number of time steps is 9, i.e., the ending time is T = 0.45. 
The material coefficient inside is ai = l.Oe -|- 06 and inside is 
ai = 1.0. The simulation results on the cutting plane (see the left plot 
in Fig. 12) is shown in Fig. 13. The relative residual error is set to 


lOe —09 as stopping criteria of the linear solvers. The iteration numbers 
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Figure 10. Simulation results with the moving object 
in the domain at different time levels t — 0.0625A:, k = 
1,9, on the cutting face. 



Figure 11. Iteration numbers (left) and CPU time 
measured in second s (right) for solving the time de¬ 
pendent heat equation with the immersed moving ob¬ 
ject in each time step: AMG preconditioned CG (solid 
lines with circle markers), AMG preconditioned GMRES 
(solid lines with star markers) in the monolithic method, 
AMG preconditioned GG (solid lines with plus mark¬ 
ers) for the Schur complement equation in the segregated 
method. 
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and the computational CPU time (in second) of the AMG precondi¬ 
tioned CG and GMRES methods in the all-at-once method, and the 
iteration numbers of the AMG preconditioned CG for the Schur com¬ 
plement equation and the computational CPU time in the segregated 
method, are shown in Fig. [T^ We observe that, in terms of iteration 
numbers, the GMRES method shows the best performance, then the 
CG method, and last the segregated method. However, regarding CPU 
time, we see that again, the CG method shows the best performance, 
then the segregated method, and last the GMRES method. 



Figure 12. Cutting plane (left), constructed growing 
sphere surfaces at time t = 0 (middle) and t = 0.45 
(right). 


6. CONGLUSION 

In this work, we develop an ALE method on the underlying ref¬ 
erence domain decomposed macro-elements consisting of tetrahedral 
and octahedral elements. That is combined with the interface cap¬ 
turing method. The numerical results demonstrate the robustness of 
this method with respect to large displacement or deformation of the 
moving interface in the model parabolic problem. We have compared 
the algebraic multigrid based all-at-once and the segregated methods 
for solving the linear system of algebraic equations arising from the 
finite element discretization. We observed that the all-at-once AMG 
preconditioned CG method shows the best performance in terms of 
CPU time. The segregated method shows comparable performance. 
Regarding the iteration numbers, the AMG preconditioned GMRES 
method shows the best performance. 
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Figure 13. Simulation results with the growing object 
in the domain at different time levels t = 0.05/c, k = 
1, 9, on the cutting face. 
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